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1. Introduction 

Harmonic and generalized harmonic (GH) coordinates have played important roles in 
general relativity theory from the very beginning. Einstein used harmonic (then called 
isothermal) coordinates in his analysis of candidate theories of gravitation (as recorded 
in his Zurich notebook of 1912) before general relativity even existed p], DeDonder 
used them to analyze the characteristic structure of general relativity in 1921 mm, 
and Fock used them to analyze gravitational waves in 1955 [3]. Harmonic coordinates 
played an important role in the proofs of the well-posedness of the Cauchy problem for 
the Einstein equations by Choquet-Bruhat in 1952 0111 and by Fischer and Marsden 
in 1972 [7]. Harmonic coordinates have also been used to obtain numerical solutions 
of Einstein’s equations by Garfinkle [8] and by Winicour and collaborators [SEMI]. 
The idea of specifying arbitrary coordinate systems using a generalization of harmonic 
coordinates was introduced by Friedrich in 1985 |T2j. And quite recently the GH 
approach to specifying coordinates played an important, perhaps seminal, role in the 
state-of-the-art numerical simulations of the final inspiral and merger of binary black- 
hole systems by Pretorius [03] jT4] using a form of the equations suggested by Gundlach, 
et al. [15], 

We think there are two important properties that have made harmonic or GH 
coordinates such an important tool throughout the history of general relativity theory. 
The first property is well known: this method of specifying the coordinates transforms 
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the principal parts of the Einstein equations into a manifestly hyperbolic form, in 
which each component of the metric is acted on by the standard second-order wave 
operator. The second property is not as widely appreciated: this method of specifying 
coordinates fundamentally transforms the constraints of the theory. This new form 
of the constraints makes it possible to modify the evolution equations in a way 
that prevents small constraint violations from growing during numerical evolutions— 
without changing the physical solutions of the system and without changing the 
fundamental hyperbolic structure of the equations. The purpose of this paper is to 
explore and understand these important properties and to extend the GH evolution 
system in a way that makes it even more useful for numerical computations. In Sec. [2] 
we review the modified form of the GH evolution system of Gundlach, et al. and 
Pretorius. We convert and extend this system in Sec. [3] into a symmetric-hyperbolic 
first-order evolution system that has constraint suppression properties comparable 
to those of the second-order system. We derive and analyze the well-posedness of 
constraint-preserving and physical boundary conditions for this new first-order system 
in Sec. [2 and in Sec. [5] we present numerical tests that demonstrate the effectiveness 
of its constraint suppression properties and the new constraint-preserving boundary 
conditions. 

2. Generalized Harmonic Evolution System 

Harmonic (sometimes called wave) coordinates are functions x a that satisfy the 
covariant scalar wave equation. These coordinates are very useful because they 
significantly simplify the second-derivative terms in the Ricci curvature tensor. To 
see this explicitly, consider a spacetime with metric tensor ip a b- 

ds 2 = ijjabdx a dx b . (1) 

(We use Latin indices from the first part of the alphabet a, b, c, ... to denote 4- 
dimensional spacetime quantities.) A coordinate x b is called harmonic if it satisfies 
the scalar wave equation, 

o = ip ab v c v c x b = -r 0 , (2) 

where V c denotes the covariant derivative compatible with ip a b, and T a = ip bc T abc is 
the trace of the standard Christoffel symbol T a t, c : 

T'abc — 2 (^b^Pac T ^c^Pab ^a^bc)- (3) 

The right side of Eq. © is just the expression for this covariant wave operator acting 
on x b in terms of partial derivatives and Christoffel symbols. 

The Ricci curvature tensor can be written as 

Rab = - ^ Cd d c dd^ab + V( a T fe ) + 1 p cd 1 p e ^ (d e l/j ca df1pdb - F ac e r bdf) , 

(4) 

in any coordinate system, where VaTf, = daR b — ij) cd T cab Td- In harmonic coordinates, 
r o = 0, so the only second-derivative term remaining in the Ricci tensor is %j> cd d c ddip a b- 
Therefore in harmonic coordinates the vacuum Einstein equations, R ab = 0, form a 
manifestly hyperbolic system [5], 

r d d c dd^ab = 2 l/j Cd 1p ef (delpcadfl/jdb - ^ace^bdf) ■ (5) 

Friedrich 12] (and independently Garfinkle 0) realized that the manifestly 
hyperbolic form of the Einstein system, Eq. m, can also be achieved for arbitrary 
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coordinates, if the choice of coordinates is fixed in a certain (but non-standard) 
way. This alternate method of specifying the choice of coordinates, which we call the 
generalized harmonic (GH) method, is implemented by assuming that the coordinates 
satisfy the inhomogeneous wave equation, 

H a {x, VO = *pabV c V c x b = -r„, (6) 


where H a (x,il>) is an arbitrary but fixed algebraic function of the coordinates x a and 
the metric ip ab (but not its derivatives). In these GH coordinates H a = — T a , so the 
vacuum Einstein equations are again manifestly hyperbolic: 

ip cd d c d d if ab = - 2 V {a H b) + 2 V' cd V' e/ (d e i/j ca d f ^ db - T ace T bdf ). (7) 

The term containing H b on the right side of Eq. 0 is a pre-specified algebraic function 
(of x a and ip a b) that operates as a source term, rather than one of the principal terms 
containing second derivatives of ip ab . The principal (i.e., second-derivative) parts of 
this GH evolution system, Eq. 0, are therefore identical to those of the harmonic 
evolution system, Eq. ©■ 

To understand the GH method of specifying coordinates more clearly, it is helpful 
to compare it to the more traditional way of specifying coordinates with the lapse 
and the shift. To do this we introduce a foliation of the spacetime by spacelike 
hypersurfaces, and adopt a coordinate system, {7, rr fc }, with the t = constant surfaces 
being the leaves of this foliation. The traditional lapse N, shift N k , and 3-dimensional 
spatial metric gij associated with this coordinate system are then defined by 

ds 2 = V 1 ab dx a dx b = — N 2 dt 2 + gij(dx l + N l dt)(dx : ’ + N^dt). (8) 


(We use Latin indices to denote 3-dimensional spatial quantities; while 

Latin indices from the first part of the alphabet a, 6, c,... will continue to denote 
4-dimensional quantities.) Expressing the GH coordinate condition, Eq. 0, in this 
3+1 language implies evolution equations for the lapse and shift: 

d t N - N k d k N = — N[H t — N l Hi + NK) , (9) 

d t N i - N k d k N l = Ngv f N(Hj + g kl T jkl ) - 8 3 n] , (10) 


where K is the trace of the extrinsic curvature. Specifying the GH gauge function 
H a (x,ip) therefore determines the time derivatives of the lapse N and shift N k , and 
hence the evolution of the gauge degrees of freedom of the system. Some gauge 
conditions (e.g., N = 1, N k = 0) may not be simple conditions on H a , just as some 
gauge conditions (e.g., H a = 0) are not simple conditions on N and N k . In this paper 
we restrict attention to the cases where H a (x, ip) is a specified algebraic function. Any 
chosen coordinates can clearly be described (ex post facto ) by an H a of this form. But 
H a may also be specified in more general ways, e.g., by giving evolution equations for 
H a I13j- We expect (but have not proven) that any coordinates can be obtained by 
specifying a priori suitable (possibly complicated) conditions on H a . 


2.1. Constraint Evolution 

Our experience in solving the Einstein equations numerically is that small constraint 
violations typically grow into large constraint violations that quickly make the 
solutions unphysical. We think it is essential therefore to understand the constraints 
and how violations of those constraints evolve with time. To this end it is helpful to 
consider the following representation of the GH system, Eq. 0: 

0 = Rab - V (a C b ), 


( 11 ) 
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where R a & is the Ricci tensor defined in Eq. ©, and C a is defined as 

C a = H a + T a . (12) 

From this perspective the condition C a = 0 serves as the constraint that ensures the 
coordinates satisfy the GH coordinate condition, Eq. 0. It is straightforward to 
verify that Eq. ED is equivalent to the GH evolution equations, Eq. 0. This form 
of the GH system, Eq. ED: is also formally equivalent to the Z4 system m ( in the 
sense that there is a one-to-one correspondence between solutions of the two systems), 
where the constraint C a plays the role of the Z4 vector field m- The systems differ 
however in the way the fields are evolved: in the Z4 system the field C a is evolved as a 
separate dynamical field, while in the GH representation C a is treated as a constraint 
which is not evolved separately. 

The evolution equation for the constraints is easily deduced from the GH evolution 
system, Eq. ED: take the divergence of the trace-reversed Eq. ED: use the contracted 
Bianclii identity X/ a R ab — ,R = 0, and exchange the order of covariant derivatives 
with the Ricci identity, yielding 

0 = \7 b \7 b C a +R ab C b . (13) 

Finally the Ricci tensor can be eliminated using Eq. ED to produce the following 
equation for the evolution of the constraints m 

0 = V b V b C a +C b V {a C b) . (14) 

This equation guarantees that the constraints C a will remain zero within the domain 
of dependence of an initial surface on which C a = dtC a = 0. Thus the GH evolution 
system is self-consistent. 

The standard Hamiltonian and momentum constraints of general relativity are 
encoded in the constraints of the GH system in an interesting way. Let t a denote the 
unit timelike normal to the t = constant surfaces of the foliation used in Eq. @. The 
standard Hamiltonian and momentum constraints are combined here into the single 
4-dimensional momentum constraint A4 a , which is given by the contraction of t a with 
the Einstein curvature tensor: 

M a = (Rab~ \^abR)t b . (15) 

Using Eq. ED for a spacetime that satisfies the GH evolution system, we see that 

t b V b C a = 2M a + ( g bc t a - t c g b a )V b C c , (16) 

where g ab = ipab + t a t b is the intrinsic metric to the t = constant hypersurfaces. 
Specifying the initial data needed to determine the evolution of the constraints, 
{C a ,d t C a }, via Eq. (THl) is equivalent therefore to specifying the more usual 
representation of the constraints, {C ai M a }, on that surface. 

2.2. Constraint Damping 

The impressive numerical simulations of binary black-hole spacetimes performed 
recently by Pretorius m m are based on a modified form of the GH evolution 
system suggested by Gundlach, et al. m ■ This modified system has the remarkable 
property that it causes constraint violations to be damped out as the system evolves. 
The modified system is given by 

0 — Dab ^(a^b) 4“ TO [^(o^ b) 2^ a b ^ ^ c ] : 


(17) 
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where t a (as before) is the future directed timelike unit normal to the t = constant 
surfaces, and 70 is a constant that determines the timescale on which the constraints 
are damped. This system can also be written more explicitly as 

il) cd d c d d ip ab = - 2X7 {a H b) + 2i/j cd ip ef (d e il) ca dfilj db - T ace T bdf ) 

+ 7o [2 6%t b) - i> ab t c ] (H c + T c ). (18) 

This system is manifestly hyperbolic since the additional constraint damping terms 
(i.e., those proportional to 70 ) do not modify the principal parts of the standard GH 
evolution system. It is also clear that the constraint-satisfying solutions of this system 
are identical to those of the standard Einstein system. 

In order to understand how this modification affects the constraints, we must 
analyze the associated constraint evolution system. This can be done by following the 
same steps that lead to Eq. m, but in this case we obtain 

0 = V b V b C a + R ab C b — 2 7o V b [ t (b C a) ], (19) 

or using Eq. m, 

0 = V c V c C Q - 2 7o V b [;£ (b C a) ] + C b V (a C b) - i 7o t a C b C b . (20) 

This constraint evolution system has the same principal part as the unmodified system, 
Eq. (1141) . Therefore the same arguments about the self-consistency of the system and 
the preservation of the constraints within the domain of dependence apply. Similarly 
the relationship between the C a constraint and the standard 4-dimensional momentum 
constraint is not changed in any essential way: setting C a = dtC a = 0 on a t =constant 
surface is still equivalent to setting C a = M a = 0 there. 

Consider the properties of the constraint evolution system for states that are very 
close to the constraint-satisfying submanifold C a = dtC a = 0 . We can ignore the terms 
in Eq. (1201) that are quadratic in C a in this case, so the constraint evolution system 
reduces to 

0 = V b V b C a ^2 7 o V b [t (b C a) ]. (21) 

Gundlach, et al. [15] have shown that all the “short wavelength” solutions to this 
constraint evolution system are damped at either the rate e~ lot or e -7ot / 2 . This 
explains how the addition of the terms proportional to 70 in the modified GH system, 
Eq. (1171) . tend to damp out small constraint violations. This also explains (in part) why 
the numerical evolutions of this system by Pretorius were so successful. A complete 
understanding of how the long wavelength constraints are damped (or not) in generic 
spacetimes would also be quite interesting, but this is not yet fully understood. 

3. New First-Order GH Evolution System 

In this section we present a new first-order representation of the modified GH evolution 
system, which will (we think) be a useful counterpart to the second-order system 
described in Sec. 12.21 above. There is an extensive mathematical literature on first- 
order evolution systems that clarifies numerous issues of great importance in numerical 
relativity, e.g., how to formulate well-posed boundary conditions [Bill HD], which 
systems form shocks GH, etc. We have also been more successful implementing first- 
order systems in our spectral evolution code. 

The principal part of each component of the modified GH system, il> cd d c d d il)ab, 
is the same as the principal part of the covariant scalar-field system. So a first-order 
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representation of the GH system can be constructed simply by adopting the methods 
used for scalar fields [22] H3]- Using this method and the usual 3+1 coordinates, 
Eq. ©, a first-order representation of the GH system can be written down, and 
indeed was written down (in essentially this form) by Alvi [21] : 

dtipab - N k d k ^ ab ~ 0, (22) 

d t n ab - N k d k n ab + Ng ki d k $ iab ~ 0, (23) 

dt^ab - N k d k $ iab + NdiU ab ~ 0, (24) 

where &i ab = diip ab and n a {, = — t c d c il> ab are new fields introduced to represent the 
first derivatives of ip a b- The notation ~ indicates that only the principal parts of the 
equations (i.e., the parts containing derivatives of the fields) are displayed. 

In the discussion that follows, it will be helpful to discuss first-order evolution 
systems like this using a more compact and more abstract notation. Thus, we let 
u a = {ipab, n a {>, &iab} denote the collection of dynamical fields; and the evolution 
system for these fields can be written as 

d t u a + A ka pd k uP = F a , (25) 

where A ka $ and F a may depend on u a but not its derivatives. We use Greek indices 
throughout this paper to label the collection of dynamical fields. The principal 
part of this system is written abstractly as d t u a + A ka pd k u& ~ 0, so Eqs. (1221) 
(1241) determine the matrix A ka p but not F a for this system. First order evolution 
systems of this form are called symmetric hyperbolic if there exists a symmetric positive 
definite matrix S a p (the symmetrizer) on the space of fields that satisfies the condition 
S a/ j,A kll p = A k p = Ap a . The mathematical literature on symmetric hyperbolic 
systems is extensive, and includes for example strong existence and uniqueness 
theorems la mum Eo]. Alvi’s representation of the GH system [24] is symmetric 
hyperbolic, as was a similar representation of the Einstein system (for the case of 
harmonic coordinates) given earlier by Fischer and Marsden [7]. 

Alvi’s first-order representation of the GH system has two serious problems: First, 
the use of the field $i ab introduces a new constraint, 

Ciab = ^iab • (26) 

which can (and does) tend to grow exponentially during numerical evolutions. Second, 
this system does not satisfy the mathematical condition (linear degeneracy) that 
prevents the formation of shocks from smooth initial data [21] . The principal part 
of the ipti component of Eq. 422D, for example, can be written as d t N l — N k d k N l ~ 0; 
and these terms have the same form as those responsible for shock formation in the 
standard hydrodynamic equations. 

We had previously developed ways to modify systems of this type to eliminate 
either of these problems [23]. However, these methods produce systems that are not 
symmetric hyperbolic when both problems are corrected simultaneously. Here we 
present new modifications that solve both problems without destroying symmetric 
hyperbolicity. We do this by adding appropriate multiples of the constraint Ci ab to 
each of the equations: ^ nN z Ci a b to Eq. (1221) . ^zN’Ciab to Eq. (1231) . and 72 NCi ab to 
Eq. (l24l) . These terms modify the principal parts of the equations: 

dti>ab - (1 + 71 )N k d k tp ab ~ 0, 

d t n ab - N k d k n ab + Ng ki d k $ iab - j 3 N k d k ^ ab ~ o, 

dt$iab - N k d k $iab + NdiUab - 72 Ndilpab ^ 0 . 


(27) 

(28) 
(29) 
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Choosing 73 = 7172 makes this new system symmetric hyperbolic for any values of 
the parameters 71 and 72 . The symmetrizer metric (which defines the energy norm) 
for this new system can be written as 

S a/ 3 du a du f3 = m ab m cd (A 2 dip a cd'ip bd + dH ac dH bd 

— 272 dip ac dH bd + g Zd d^iacd^jbd) , (30) 

where m ab is any positive definite metric (e.g., m ab = g ab + t a t b or even m ab = S ab ) 
and A is a constant with dimension length -1 . This symmetrizer is positive definite so 
long as A 2 > 

The eigenvectors of the characteristic matrix, n k A ka p (where n k is the outward 
directed unit normal to the boundary of the computational domain), play an important 
role in setting boundary conditions for first-order evolution systems. Let e“^ denote 
the left eigenvectors with eigenvalues v^, defined by 

e 0 ‘ ll n k A k »p = v {&) e a p. (31) 

We use indices with hats (e.g., a) to label the characteristic eigenvectors and 
eigenvalues, and a is not summed over in Eq. m■ The eigenvalues v^) are also 
called the characteristic speeds. The characteristic matrices of symmetric hyperbolic 
systems have complete sets of eigenvectors, so the matrix e a p is invertible in this case. 
The characteristic fields, u a , are defined as the projections of the dynamical fields onto 
the characteristic eigenvectors: u a = e a gu^. Boundary conditions must be imposed 
on each incoming characteristic held, i.e., each u a with negative characteristic speed, 
V(a) < 0 US EH 20]. The characteristic fields for the new GH evolution system, 
Eqs. m HMD, are given by 


Uab 

— 0 abt 

(32) 

Kb 

n a6 i Tl ^iab T2 

(33) 

^iab 

= Pi k ®kab, 

(34) 


where Pi k = 5 k — ntn k . The characteristic Helds u° ab have coordinate characteristic 
speed —(1 + r yi)n k N k , the Helds m 1 ^ have speed —n k N k ±N, and the fields u 2 ab have 
speed —n k N k . 

The complete equations for our new first-order representation of the GH evolution 
system (including all the non-principal parts) are 

dt^ab - (1 + 7i )N k d k ip ab = -NIL ab - 7iiV%a6, (35) 

d t n ab - N k d k n ab + Ng ki d k $iab - 7172 N k d k ^ ab 

= 2 N^ cd (g l ^ lca ^ jdb - n ca n db - i) ef T ace T bdf ) 

- 2NV( a H b ) — \Nt c t d n cd R ab — Nt c H ci g lj $j ab 
+ iV 7o [2 S c (a t b) - 1/J ab t c ] (H c + r c ) - 7172^$^, (36) 

dt^iab - A T k d k $iab + Ndili ab - N^di^ab 

= \Nt c t d <b icd n ab + Ngi k t c <f> ijc <t>kab - N l2 <b lab . (37) 

The terms on the right sides of Eqs. (l3lH) (Ef7l) are algebraic functions of the dynamical 
fields. The connection terms T cab appearing on the right side of Eq. (ED are computed 
using Eqs. ©, where it is understood that the partial derivatives are to determined 
from the dynamical fields by 

dtipab = ~ NH ab + N l $ iab , 

di'&ab = *&iab- 


(38) 

(39) 
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Choosing the parameter 72 > 0 in this new system causes the constraint Ci ab to 
be exponentially suppressed [23], because the modified Eq. (1371) implies an evolution 
equation for C iab having the form, d t C iab - N k d k C iab ~ - 7 2 NC iab . Choosing 7 ! = -1 
makes the system Eqs. (l35l) (1371) linearly degenerate, which implies that shocks do not 
form from smooth initial data EH- And choosing the parameter 70 > 0 causes the 
constraint C a to be exponentially suppressed, as discussed in Sec. 12.21 

4. Boundary Conditions 

The modifications of the GH evolution system discussed in Secs. l2.2l andl3lare designed 
to damp out small constraint violations that may arise from inexact initial data, 
numerical errors, etc. These modifications will do nothing, however, to prevent the 
influx of constraint violations through the boundaries of the computational domain. 
Constraint-preserving boundary conditions are needed to prevent this G3 GDI E3 E3 
[25]. Such boundary conditions can be formulated once the propagation equations for 
the constraints are understood. So we derive a first-order system of evolution equations 
for the constraints in Sec. mi 1 use them to derive constraint-preserving boundary 
conditions in Sec. mi present boundary conditions for the physical gravitational-wave 
degrees of freedom of the system in Sec. 14.31 and finally analyze the well-posedness of 
the combined set of new boundary conditions in Sec. 14.41 

4-1. First-Order Constraint Evolution System 

The primary constraint of the GH system is the gauge constraint, C a , which we re-write 
here in terms of the first-order dynamical fields: 

Ca = H a + g ij <f> ija + t b n ba - yi^ibc - \t a ^ bc n 6c . (40) 

This expression differs from Eq. (fT21) only by multiples of the constraint Ci ab . In the 
following we use this definition, Eq. <@QD, rather than Eq. m, because it simplifies 
the form of the constraint evolution system. The evolution equation for C a , Eq. m, 
is second order. Thus, we must define new constraint helds that represent the first 
derivatives of C a in order to reduce the constraint evolution system to first-order form. 
Thus we define new constraint helds T a and Ci a that satisfy 

T a « t c d c C a = N-\d t C a - N%C a ), (41) 

C ia « d,C a , (42) 

up to terms proportional to the constraints C a and Ci ab . The following definitions of 
T„ and Cia accomplish this in a way that keeps the form of the constraint evolution 
system as simple as possible: 

Ta = - g ij diH ja - gWd&jba + \t a ^ bc g i3 d^ jbc 

+ t a g 13 d,Hj + gi* i j b gi k $ kcd il> bd f - \g^ ijb g jk ^kcdr d t h 

- 9 l at b diH b + g v $icd$jbai’ bc t d - lt a g IJ g mn $ imc $ n j d ip cd 

- y a g ii *icd* j beil>'*V U + \t a H cd H be ilj cb ^ de - g ij HiU ja 

- t b g ij n bi n ja - \g^ lcd t c t d u be ^ be + ±t a n cd n be -iij ce t d t b 
+ g i a $i cd n be t c t b ^ de - g i i$ iba t b U je t e - y ij $ icd t c t d n ja 

- g ii H& jba t b + g l a ^iedH b ^ bc t d + 72 (g ld C lda - 

+ \t a H cd r d H b t b - t a g i3 ^i jc H d r d + k t a9 ii H i $ jcd r d , (43) 
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c ia = g jk dj$ika - \gir%*icd + t b diU ba - \t a ^ cd diH cd 

+ d,H a + y a <s> 0Cd <s> ief r e Vf + y k ^jcd^ike^ cd t e t a 

- g jk g mn $ jm a$ikn + ^icdliheta (V^V^ + ^t°t d ) 

- $ lcd n ba t c (ip bd + \t b t d ) + i 72 (t a r d - m d ) C lcd . (44) 

The remaining constraints needed to complete the GH constraint evolution system are 
Ciab defined in Eq. (EiD, and the closely related Cij ab , defined by 

Cijab — = 2 d[jCj] ab . (45) 

The complete collection of constraints for the GH evolution system is therefore the 
set c A = {C a ,E a ,C ia ,C iab ,Cijab} defined in Eqs. (00]), J43]), (|44j) . ([26]) : and ([45]). (We 
use upper case Latin indices to label the constraints.) The constraints c A depend on 
the dynamical fields u a = {ip a b, n o b, d>i a t} and their spatial derivatives dku a . Thus 
the evolution of the constraint fields c A is completely determined by the evolution 
of the dynamical fields through Eqs. m iii)- We have evaluated these constraint 
evolution equations and have verified that they can be written in the abstract form 

d t c A + A kA B (u)dkC B = F A B (u,du) c B , (46) 

where A kA B and F A B may depend on the dynamical fields u a and their spatial 
derivatives dkU a . Thus the constraint evolution system closes: the time derivatives 
of the constraints vanish initially when the constraints themselves vanish at an initial 
time. The principal part of the first-order constraint evolution system turns out to 
be remarkably simple (given the complexity of the expressions for the constraints 
themselves): 


dtC a 

— o, 

(47) 

d t T a 

~ IPdiTa + Ng^diCja, 

(48) 

dtCia 

- N^djCia + Nd.Ta, 

(49) 

&tCiab 

~ (1 + ^)N k d k C iab , 

(50) 

^tCijab 

~ N k d k C ijab . 

(51) 


This constraint evolution system is symmetric hyperbolic with symmetrizer 

S AB dc A dc B = m ab [dT a dTb + g ij (. dCiadC jb + g kl m cd dC ikac dC jlbd ) 

+ A 2 (dC a dC b + g lj m cd dCiacdCjbd) , (52) 


where A 2 is a positive constant and m ab is an arbitrary positive definite metric. The 
constraint energy for this system is defined as 


Co 


Sa B c a c b yfgd^x. 


(53) 


Since the constraint evolution system is hyperbolic, it follows (at the continuum level) 
that the constraints will remain satisfied within the domain of dependence of the initial 
data, if they are satisfied initially. 

We have analyzed the solutions to this constraint evolution system for the case 
of small constraint violations of solutions near flat space. We find that all of the 
short-wavelength constraint violations are damped at the rate e~ lot , e -70 */ 2 , or 
e -72t . So choosing 7o > 0 and 72 > 0 is sufficient to guarantee that all of these 
constraints are suppressed. This new first-order GH system therefore has the same 
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constraint suppression properties as the second-order system of Gundlach, et al. |15] 
and Pretorius m- 

The constraint evolution system, Eq. (1461) , is symmetric hyperbolic and it will be 
useful to determine the characteristic constraint fields. Thus, we evaluate the matrix 
of left eigenvectors of the constraint evolution system e A b and their corresponding 
eigenvalues v,^ (or characteristic speeds). The characteristic constraint fields are 
defined (in analogy with the principal evolution system) as the projections of the 
constraint fields onto these eigenvectors: c A = e A BC B . The resulting characteristic 
fields for this constraint evolution system are 


cr 

— 'Fa “F ^ Ckai 

(54) 

Co 

= C a , 

(55) 

c 2 

^ia 

= P k iC ka , 

(56) 

3 

0 iab 

C'iab'i 

(57) 

J 

^ijab 

— Cijab- 

(58) 


The characteristic constraint fields c„ have coordinate characteristic speeds — mN 1 ± 
N, the fields c\ have speed 0, the fields cf a and c A jah have speed —mN 1 , and the fields 
c iab lmve speed -(1 + 7 i )niN l . 


4-2. Constraint-Preserving Boundary Conditions 


Boundary conditions must be imposed on all the incoming characteristic fields u a , 
i.e., all those with v^) < 0 on a particular boundary. Thus, boundary conditions 

will typically be needed for the characteristic field uz7, and (depending on the value 
of the parameter 71 and the orientation of the shift N k at the boundary) may also 
be needed for u^ b and/or u~ ab . Some of these boundary conditions must be set by 
physical considerations, i.e., by specifying what physical gravitational waves enter the 
computational domain. Some of the boundary conditions can be used, however, to 
prevent the influx of constraint violations. This can be done by specifying the incoming 
u a at the boundary in a way that ensures the incoming characteristic constraint fields 
c A also vanish there. The incoming constraint fields for this system include c° - , and 
perhaps c/ ab and/or cf kab depending on 71 and N k at the boundary. We find that 
these incoming c A are related to the incoming u a by the following expressions: 


C°- « V2 \k {c ip d) c 


i n 3 


n c, 


iab 


d±u 


^ ^ikab ~ 


ab") 

2 



d±u 


1 - 
cd ’ 


(59) 

(60) 
(61) 


Here the notation d±u a denotes the characteristic projection of the normal derivatives 
of u a (i.e., d±u a = epii k dkvP ), and ss implies that algebraic terms and terms involving 
tangential derivatives of the fields (i.e., P k idkU a ) have not been displayed. The inward 
directed null vector k c used here is defined as k c = ( t c — n c ) / y/2. The idea is to set 
the left sides of Eqs. (ESI) ED to zero to get Neumann-like boundary conditions for 
the indicated components of d±u a . By imposing these conditions on d±u a , we ensure 
that these incoming components of c A vanish. 
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We have found that a convenient way to impose boundary conditions of this type 
is to set the incoming projections of the time derivatives of u a , d t u a = e“ pd t uP , in 
the following way: 

d t u a = D t u a + V( & ) (d±u a - d±u a \ gc ). (62) 

In this expression the terms D t u a represent the projections of the right sides of the 
evolution system, Eqs. © (EZD; so the equations at non-boundary points would 
simply be dtu a = D t u a . The term d±u a \ BC is the value to which d±u a is to 
be fixed on the boundary. This form of the boundary condition replaces all of the 
d±u a that appears in D t u a with d±u a \ BC ,. Applying this method to the constraint¬ 
preserving boundary conditions in Eqs. (|59| HMD, we obtain the following rather 
simple conditions 

dtu° ab = D t u° ab - (1 + 7 i )njN J n k cl ab , (63) 

d t ul b = [\P ab P cd -2l {a P b) H^ + l a l b k c k d ]D t u\ d 

+ V2(N + n j Ni)[l (a P b f-±P ab l c -±l a l b k c ]c°-, (64) 

dtulab = D tu 2 kab - niN l n l P J k c A ijab . (65) 

The quantity P ab in these expressions is the projection tensor, P ab = il)ab + t a tb — n a n b , 
and the outgoing null vector l a is defined by l a = ( t a + n a )/\/2. 

4-3. Physical Boundary Conditions 

The constraint-preserving boundary conditions presented in Eqs. (P| p| restrict 
only four degrees of freedom of u A b . Two of the remaining degrees of freedom represent 
the physical gravitational waves, and the final four represent gauge freedom. We 
choose to characterize and control this gravitational wave freedom in terms of the 
incoming parts of the Weyl curvature. The propagating components of the Weyl 
tensor can be written as 

*4 = ( Pa C Pb d - \PabP cd ) ( t e T n e ) {t f T n f )c cedf . (66) 

We showed in Ref. [28] that these components of the Weyl tensor are the incoming 
and outgoing (respectively) characteristic fields of the curvature evolution system that 
follows from the Bianchi identities. The w± b are proportional to the Newman-Penrose 
curvature spinor components 4/4 (outgoing) and 4/o (ingoing) respectively. We also 
note that the spatial components of w ^ are equal to the components of the Weyl tensor 
characteristic fields 2LR defined in our paper on constraint-preserving boundary 
conditions for the KST system [28|. The expression for the Weyl tensor in terms 
of our first order variables is unique only up to terms proportional to constraints; it 
is possible to choose these constraint terms so that the depend on the normal 
derivatives of ti" in the following way: 

<4 « ( Pa C Pb d \PabP Cd ) (d±l# + 72d ± l&) . (67) 

Thus a physical boundary condition can be placed on the relevant components of u^ b 
using the method of Eq. m by setting 

dtu\ b = ( P a C Pb d - \PabP Cd ) X 

[Dtulj ~(N + njN J )(w~ d - 72 n*c| cd )]. (68) 


Generalized Harmonic Evolution System 


12 


We can also inject incoming physical gravitational waves with a predetermined 
waveform h a t(t , x) through the boundary of the computational domain by setting 

dtuli = ( P a c P b d - \P ab P cd )h cd (t , x). (69) 

The case h ab = 0 corresponds to an isolated system with no incoming gravitational 
waves. 

More generally we can combine the constraint-preserving, physical no-incoming 
radiation, and the injected gravitational wave boundary conditions by setting dtu *7 
equal to the sum of the right sides of Eqs. (El, ©, and EH), and setting the time 
derivatives of the other incoming fields according to Eqs. E3) and (fTTSl) . Note that 
this set of combined boundary conditions holds the pure gauge components of vrZ 
constant in time; other boundary conditions on the gauge degrees of freedom are of 
course possible but are not considered here. 

4-4- Well-posedness 

The well-posedness of the initial-boundary value problem can be analyzed using the 
Fourier-Laplace technique [30,. We have applied this method to the GH system with 
the combined set of boundary conditions presented here: we treat the case of high- 
frequency perturbations of flat spacetime in a slicing with flat spatial metric, unit 
lapse, and a constant shift that is tangent to the boundary. Applying the Fourier- 
Laplace technique to this case yields a necessary (but not sufficient) condition for well- 
posedness, the so-called determinant condition [30| : failure to satisfy this condition 
would mean the system admits exponentially growing solutions with arbitrarily large 
growth rates. We have verified that this determinant condition is satisfied for the GH 
system using the combined set of boundary conditions presented here. 

5. Numerical Results 

In this section we describe several numerical tests of the new first-order GH evolution 
system. First we test the effectiveness of the two constraint damping terms included 
in Eqs. (f35l) (l37l) by evolving Schwarzschild initial data (in Kerr-Schild coordinates). 
These tests are performed on a computational domain consisting of a spherical shell 
that extends from r m i n = 1.8 M (just inside the event horizon) to r max = 11.8 M, 
where M is the mass of the black hole. In these evolutions we “freeze” the values 
of the incoming characteristic fields to their initial values by setting dtu a = 0 on the 
boundaries for all incoming fields (i.e., all u a with < 0). We performed these 
numerical evolutions using spectral methods as described for example in Ref. [28] for 
a range of numerical resolutions specified by the parameters N r (the highest radial 
spectral basis function) and L max (the highest spherical-harmonic basis function). 
Figure |T] illustrates the results of these tests for several values of the constraint 
damping parameters 70 and 72 . These tests show that without constraint damping the 
extended GH evolution system is extremely unstable, but with constraint damping the 
evolutions of the Schwarzschild spacetime are completely stable up to t = 10, 000 M 
(and forever, we presume). These tests also illustrate that both the 70 and the 72 
constraint damping terms are essential for stable evolutions. 

Constraint violations in Fig. [T] (and in the rest of this paper) are measured with the 
constraint energy £ c defined in Eq. (l53l) . Since £ c is not dimensionless, its magnitude 
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t/M t/M 


Figure 1 . Evolution of constraint violations for Schwarzschild initial data. Left 
figure shows evolutions using various values of the constraint damping parameters 
70 and 72 using numerical resolution {N r ,Lmax} = {13,7}. Right figure shows 
the long timescale evolution of the same data for three different numerical 
resolutions. 


has no absolute meaning. We construct an appropriate scale with which to compare 
S c by evaluating the L 2 norm of the spatial gradients of the dynamical fields, 

\\du\\ 2 = J g lJ m ab m cd (^A 2 di'ip ac d j 4> bd +d l U ac djU bd 

+ g kl di$kacdj$ibd)Vgd 3 x. (70) 

The dimensionless constraint norm |\C 11 shown in these figures is defined as 


l|C|| 


w 


(71) 


which is a meaningful measure of the relative size of constraint violations in a particular 
solution. In the figures shown here we evaluate ||C|| with m ab = 5 ab and the 
dimensional constant A = 1/M. 

Our second numerical test evolves the somewhat more challenging initial data for 
a Kerr black hole (in Kerr-Schild coordinates) on a computational domain consisting 
of a spherical shell that extends from ?' m i n = 1.8 M (just inside the event horizon) 
to 7'max = 21.8 M. We use two subdomains, each having numerical resolution 
{N r , L m ax}, to cover this region. The spin of the Kerr spacetime used here is 
a = (0.1,0.2,0.3)M, where the magnitude of this vector determines the Kerr spin 
parameter a = |a| « 0.374M, and the direction determines the orientation of the 
Kerr rotation axis relative to the quasi-Cartesian coordinate system used in our code. 
For this test we use the combined set of physical and constraint-preserving boundary 
conditions discussed at the end of Sec. 14.31 Figure [2] shows that numerical evolutions 
of this Kerr spacetime are stable and numerically convergent to t = 10, 000M (and 
forever, we presume) using a range of numerical resolutions. 

Our third numerical test is designed to demonstrate the effectiveness of our 
new constraint-preserving boundary conditions. This test consists of evolving a 
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Figure 2. Evolution of constraint violations for Kerr initial data with spin 
parameter a = (0.1, 0.2,0.3) for several numerical resolutions. 


black-hole spacetime perturbed by an incoming gravitational wave pulse. We start 
with Schwarzschild initial data, and perturb it via the incoming gravitational wave 
boundary condition described in Eq. (|69l) with h a b = f{t)(x a x b + y a y b — 2 z a z b ) where 
x a , y a , and z a are the components of the coordinate basis vectors, x a d a = d x , etc. 
For these evolutions we use an incoming gravitational wave pulse whose time profile is 
fit) = Ae~^ t ~ tp ' )2 Z 1 " 2 with A = 10~ 3 , t p = 60 M, and w = 10 M. This test is performed 
on the same computational domain described above for the second numerical test. 
Figure [3] illustrates the results of these tests using two types of boundary conditions: 
frozen-incoming-field (i.e., dtu a = 0 for < 0) boundary conditions (solid curves) 
and the new combined set of constraint-preserving and physical boundary conditions 
discussed at the end of Sec. 14.31 (dashed curves). The graph on the left in Fig. [3] 
shows that constraint violations converge toward zero as the numerical resolution is 
increased when the new boundary conditions are used, but not when frozen-incoming- 
field boundary conditions are used. 

The graph on the right in Fig. [3] shows the outgoing physical gravitational wave 
flux (measured on the outer boundary of the computational domain) computed using 
frozen-incoming-fields (sold curve) and the new constraint-preserving and physical 
(dashed curve) boundary conditions. These evolutions were computed with numerical 
resolution {N r ,L m ax } = {21,19}. We measure the outgoing gravitational wave flux 
with the quantity (i?*!^), which is the Weyl curvature component *F 4 averaged over 
the outer boundary of our computational domain: 

47r(i ?^ 4 ) 2 = J |* 4 | 2 d 2 V. (72) 

Here 47 ri ? 2 is the proper surface area of the boundary, and d 2 V represents the 
proper area element on this boundary. Since 4 / 4 falls off like 1 /R, this quantity 
should be independent of R (asymptotically). The dashed curve on the right in 
Fig. [3] clearly shows quasi-normal oscillations with frequency coM = 0.376 — 0.089 i 
(determined by a numerical fit to these data). This is in good agreement with 
the frequency of the most slowly damped quasi-normal mode of the black hole: 
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Figure 3. Evolution of Schwarzschild initial data perturbed by a gravitational 
wave pulse with amplitude 10 —3 . Left figure depicts constraint violations at 
various numerical resolutions, and the right figure shows ^4 averaged over the 
outer boundary of the computational domain at a single numerical resolution. 
Solid curves use freezing boundary conditions and dashed curves use constraint¬ 
preserving and physical boundary conditions. 


u>M = 0.37367 — 0.08896 i [31]. It is interesting to note that the solid curve—using 
frozen-incoming-fields boundary conditions—gives qualitatively incorrect results for 
the physical gravitational waveform, even though the level of constraint violations is 
fairly small numerically in this case. This is not surprising because the magnitude of 
constraint violations in this case is comparable to the size of the injected gravitational 
wave pulse. 
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